en 


Open Access 


Publishing in European Networks 


Published version 


ERC-OAPEN-2015, Deposit cover page 


Open Access Book Chapter 


Chapter title Exploring the Components of the Universe Through 
Higher-Order Weak Lensing Statistics 

Title in English 

Author(s) Adrienne Leonard, Jean-Luc Starck, 
Sandrine Pires and Francois-Xavier Dupé 

DOI 10.5772/51871 

Copyright holder _| The Authors 

License CC BY 3.0 license 

Language English 

Original publication 

Book title Open Questions in Cosmology 

Title in English 

Publisher Intech 

Publication date 2012 

ISBN 9789535108801 

DOI 

Editor(s) Gonzalo J. Olmo 

ERC Grant information 

Project title EC | SPARSEASTRO 

Project acronym SparseAstro 

Grant number 228261 

ERC principal Jean-Luc Starck 

investigator 


under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 683680). 


OAPEN Foundation received funding from the European Research Council (ERC) 


Chapter 2 


Exploring the Components of the Universe Through 
Higher-Order Weak Lensing Statistics 


Adrienne Leonard, Jean-Luc Starck, 
Sandrine Pires and François-Xavier Dupé 


Additional information is available at the end of the chapter 


http://dx.doi.org/10.5772/51871 


1. Introduction 


Our current cosmological model, backed by a large body of evidence from a variety of 
different cosmological probes (for example, see [1, 2]), describes a Universe comprised of 
around 5% normal baryonic matter, 22% cold dark matter and 73% dark energy. While many 
cosmologists accept this so-called concordance cosmology - the ACDM cosmological model 
— as accurate, very little is known about the nature and properties of these dark components 
of the Universe. 


Studies of the cosmic microwave background (CMB), combined with other observational 
evidence of big bang nucleosynthesis indicate that dark matter is non-baryonic. This 
supports measurements on galaxy and cluster scales, which found evidence of a large 
proportion of dark matter. This dark matter appears to be cold and collisionless, apparent 
only through its gravitational effects. 


While dark matter is largely responsible for the growth of the largest structures in the 
Universe, dark energy — dominant at late times — appears to have a negative pressure, and to 
be responsible for an accelerated expansion of the Universe [3]. It is usually parameterised 
by its equation of state parameter w = p/p, where p is the pressure associated with the 
dark energy and p is its energy density. An equation of state parameter of w = —1 would 
indicate a cosmological constant, consistent with general relativistic theory, but deviations 
from this value would suggest a rather more exotic dark energy, and might perhaps imply a 
modification to our current theory of gravity. 


Gravitational lensing — the deflection of light rays by massive objects due to gravitational 
effects, which gives rise to distortions in images of galaxies — is an ideal probe of the dark 
universe, as it can probe the evolution of the dark matter power spectrum in an unbiased 
way, and offers complementary constraints to those obtained from the CMB and other probes 
of large-scale structure. 
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Moreover, gravitational lensing probes cosmological perturbations on smaller angular scales 
than CMB studies, and is thus sensitive to the non-Gaussianity induced by the late-time 
non-linear evolution of structures such as clusters of galaxies, as well as any primordial 
non-Gaussianity arising, for example, due to inflation very early in the evolution of the 
Universe. 


Traditionally, constraints from gravitational lensing are obtained by considering two-point 
statistics of the lensing shear field, which encodes the small elliptical distortions applied to 
the images of galaxies as a result of the gravitational potential field of structures along the 
light’s path. Such two-point statistics are only weakly sensitive to the dark energy density 
parameter Q in and depend on a degenerate combination of the amplitude of the matter 
power spectrum og and the matter density parameter Oym. In addition, two-point statistics 
probe only the Gaussian part of the shear field, therefore in considering such statistics alone, 
information about nonlinear structure evolution and primordial non-Gaussianity is lost. 


Similarly, attempts to reconstruct a map of the two-dimensional projected surface mass 
density (the convergence x) and three dimensional density field have often involved the use 
of Gaussian priors to constrain the reconstruction (for example, see [4]), thus again having 
limited application in studies of non-Gaussianity. 


In this chapter, we present a review of recently developed gravitational lensing techniques 
that go beyond the standard two-point statistics, both in the arena of map-making in two 
and three dimensions and that of higher-order statistics of the shear field or, equivalently, 
the convergence field x. 


The methods we present are all based on the concept of sparse recovery, which has been 
found to be a very powerful tool in signal processing [5, 6]. Such methods are based on the 
assumption that a given image or observation can be represented sparsely in an appropriate 
basis (such as a Fourier or wavelet basis). Sparse priors using a wavelet basis have been 
used in many areas of signal processing in astronomy; of particular interest in this chapter 
will be the areas of denoising and in the reconstruction of the 3D density field from lensing 
measurements, but other applications include deconvolution and inpainting. 


There is much information to be gained by considering non-Gaussian statistics and nonlinear 
signal processing methods, and this is an exciting and active area of research. With new 
surveys such as the Euclid satellite [7] coming online within the next decade, a wealth 
of high-quality data will soon be available. These new techniques will therefore prove 
invaluable in constraining the cosmological model, and allow us to better understand the 
nature of the primary constituents of the Universe. 


2. Weak gravitational lensing 


We begin with a brief overview of weak lensing theory. For a more complete description and 
discussion of the subject, see [8, 9] 


2.1. Weak lensing theory 
The basic idea underlying the theory of gravitational lensing is that massive objects distort 


the spacetime around them, and thus bend the path of light in their vicinity. A light ray 


1 The energy density parameters are defined by Q; = p;/¢, where p; is the critical density of the Universe 


Exploring the Components of the Universe Through Higher-Order Weak Lensing Statistics 
http://dx.doi.org/10.5772/51871 


originating at angle G is deflected such that it appears to the observer to originate at an 
angular position 8, where 


B=0-a(d) =0-VYy(8), (1) 
where y is the two-dimensional deflection potential associated with the lens and « is the 


deflection angle. 


For extended sources, photons from different angular positions in the source plane are 
deflected differently, giving rise to a distortion in the observed galaxy image, which is 
described — to first order — by the Jacobian of the lens equation (1): 


_ 08 _ ep _(1-r-n -nv 
A 00 $i, 90;90; — y2 1—xK+ Yı É 2) 


where kK = 5V7 is the convergence, or dimensionless surface density, and y = y1 + iy2 = 
|y\e2"? is the complex shear, which gives rise to an anisotropic elliptical distortion in the 
lensed image. The components of the shear are related to the potential via: 


1 
n= 5 (%- 3p), 2 = Ardy, ©) 


where 0; = 0/06;, and the shear is related to the convergence through the relation: 


(0) = a 26'D(0 — 0')x(0") , (4) 
where the convolution kernel D is given by 


62 — 02 — 210,05 1 
D(0) = 2 1 ae : 
(9) || (0; — i0)? ©) 


The convergence, in turn, can be related to the 3D density contrast ô(r) = p(r)/p — 1 by 


3HiOm fP dy! TK Ww) Fx (w — w’) d[fx(w')@, w') 


OW) Ly fa) Co 


(6) 


where p is the mean density of the Universe, Ho is the Hubble parameter, Om is the matter 
density parameter, c is the speed of light, a(w) is the scale parameter evaluated at comoving 
distance w, and 


K"? sin(K!/2w), K>0 
fx(w) = 4 4, K=0, (7) 
(—K)~\/?sinh({—K]!/2w) K <0 


gives the comoving angular diameter distance as a function of the comoving distance and 
the curvature, K, of the Universe. 
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2.2. Two-point statistics 


The typical gravitational shear applied to a galaxy as a result of lensing by large-scale 
structure in the Universe — so-called cosmic shear — is of order |y| ~ 0.01. However, galaxies 
are intrinsically elliptical in shape, with a typical ellipticity of order |e| ~ 0.2 — 0.3, therefore 
the gravitational lensing effect can only be measured statistically; under the assumption 
that galaxy shapes are intrinsically uncorrelated, the mean ellipticity computed over a large 
number of sources will yield the gravitational shear: (e) ~ y. 


The most common method for constraining cosmological parameters in weak lensing studies 
is to use two-point statistics of the shear field. The power spectrum of the shear or 
convergence, P,(¢) = P}(£), can be related directly to the 3D matter power spectrum of 
density fluctuations 6 by: 


ni PR aea) o 
where 
wow) = [ao EE) oe) |] Q 


is a weighting function integrated over the probability distribution of sources p(z) in the 
sample as a function of redshift z. 


When working with real data, it is often more convenient to consider statistics computed 
in real space, namely the shear correlation functions, which are defined in terms of the 
convergence power spectrum as [10, 11]: 


E0) = (re) End = [°F Tall), (10) 


where Jn is an n-th order Bessel function of the first kind, and y+ and yx are the tangential 
and cross shear components, respectively, which are defined relative to the vector connecting 
the two galaxies separated by angular distance 0. 


2.3. Cosmological constraints 


The measured two-point statistics described above can, by virtue of their relationship to the 
power spectrum of the underlying matter density fluctuations, be used to place constraints 
on cosmological parameters. These constraints may be improved if information about the 
distances to source galaxies is used to bin the galaxies into redshift slices, and compute the 
correlation functions tomographically. 


Figure 1 shows cosmological parameter constraints obtained using two state-of-the-art 
telescopes: the Canada-France-Hawaii Telescope (CFHT), and the Hubble Space Telescope 
(specifically, the Cosmos survey). While the CFHT results were obtained using a 2D analysis 
[12], the Cosmos results show constraints obtained by both a 2D and a 3D analysis of 
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the shear data [13]. In both cases, a large degeneracy is seen between the matter density 
parameter Om and og, the amplitude of the matter power spectrum. The CFHT results 
also show constraints on the dark energy equation of state parameter w. Here, again, a 
degeneracy between parameters is seen. 
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(a) Joint constraints on Om and og (b) Joint constraints on Om and w (c) Joint constraints on Om and gg from 
from a 2D analysis of the combined from a 2D analysis of the CFHT deep a 2D and 3D analysis of the Cosmos 
CFHT wide and deep surveys [12]. survey [12]. survey data [13]. 


Figure 1. Cosmological parameter constraints from recent weak lensing surveys. Contours are plotted at the 68% (cyan), 95% 
(blue) and 99.9% (green) confidence levels for the CFHT data, and at the 68% (inner) and 95% (outer) confidence levels for the 
Cosmos data. 


While it is clear from the figures that two-point shear statistics hold a wealth of information, 
it is also evident that these statistics depend on degenerate combinations of the cosmological 
model parameters. In order to break these degeneracies, and thus more tightly constrain 
our cosmological model, we must therefore consider higher-order statistics of the shear 
or convergence. Such higher-order statistics will enable us to capture information about 
the non-Gaussian part of the lensing spectrum due to a combination of primordial 
non-Gaussianity and late-time nonlinear evolution of structures, which two-point statistics 
are unable to probe. 


3. Scalar fields associated with the shear 


The shear y is a spin-2 field, therefore while the two-point correlations of the shear field 
can be reduced to a scalar quantity for parity reasons, this is not the case for higher-order 
moments of the shear field [14-16]. For this reason, it is useful to compute a scalar quantity 
from the shear field before computing higher-order statistics. 


3.1. The aperture mass statistic 


For this purpose, the aperture mass statistic Map [17, 18] is widely-used. The Map statistic 
is defined as the convolution of the convergence x with a radially symmetrical filter function 
U(|#|) of width 8: 


Map(@) = [ve K(B)U(|9|) . (11) 
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By considering the relationship between the shear, y, and the convergence, one can 
reformulate equation (11) in terms of the measured shear as 


Map(6) = | 28 tOO), a2) 


where y:(®) is the tangential component of the shear at position ® relative to the centre of 
the aperture, Q(|?|) is a second radially-symmetric function, related to U(|#|) by: 


Damn fit? 
Q(8) = zl PUJAS UCO), (13) 
and U(#) is required to be compensated, i.e. 
Oeut 
j 0 U(8) dd =0, (14) 
JO 


with cy; often taken to be the radius of the aperture, 0. Furthermore, Q(@) and U(#) are 
required to go to zero smoothly at Oy. It is also preferable that the power spectrum of U(#) 
is local in the frequency domain, and shows no oscillatory behaviour. This ensures that the 
filter function acts as a band-pass filter, allowing detection of structures at the scale of interest 
only. 


Several authors [19, 20] have advocated a filter function of the form: 


A bv be? 
u) = z ( x ) oP ( 7) ; 


A b8? bv 
Q(0) = @2 62 exp ( 02 ) , (15) 


where various choices for the constant b and the overall normalisation A have been used in 
the literature [20-22]. Of specific interest is the form used in [20] and [22], where b = 1/2 is 
chosen.? This form is considered to be optimal for higher-order weak lensing statistics, such 
as the skewness of the aperture mass statistic [19, 22]. 


Note that the Q(¥) filter function described in equation (15) shows a peak at 8 = /20, and 
tends to zero as Ŷ — ov. In practice, any algorithm used to generate an aperture mass map 
will need to truncate this filter function at some finite radius. This will involve a trade-off 
between accuracy and algorithm speed, and truncation may affect the effective width of the 
filter. 


One advantage of this method is that the aperture mass statistic can be computed directly 
from the shear catalogue. Moreover, as the filter acts as a band-pass filter, it is possible to 
boost the signal relative to the noise by considering a filter with a scale substantially larger 


? This form is analogous to the Mexican Hat wavelet function 
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than the typical scale of the noise (which is typically dominant on pixel scales, depending 
on the binning chosen for the shear data). Indeed, optimal signal-to-noise is obtained when 
the filter is chosen such that its angular scale matches that of the structures we aim to detect, 
and its shape matches closely the expected profile of these structures. 


3.2. The wavelet transform 


The wavelet transform is a multiscale transform, where the wavelet coefficients of an image 
are computed at each position in the image at various different scales simultaneously. In one 
dimension, the wavelet coefficient of a function f(x), evaluated at position b and scale a is 
defined as [23, 24]: 


wab) =e f tow (EF) ax, (16) 


where (x) is the analysing wavelet. The analysis is analogous in 2 dimensions, with 
p(x y) = P(x) ply). 


By definition, wavelets are compensated functions; ie. the wavelet function p(x) is 
constrained such that Jip: Y(x)dx = 0 and hence, by extension 


i (x, y)dx dy = 0. (17) 


According to the definition in equation (16), the continuous wavelet transform of an image is 
therefore nothing more than the convolution of that image with compensated filter functions 
of various characteristic scales. If the image f(x,y) is taken to be the convergence (x,y), 
then for an appropriate choice of (radially-symmetric, local) wavelet, the wavelet transform 
is formally identical to the aperture mass statistic at the corresponding scales, the only 
difference being the choice of filter functions. 


In practice in application, we use the starlet transform algorithm [6, 23-25], which 
simultaneously computes the wavelet transform on dyadic scales corresponding to 2j pixels. 
This algorithm decomposes the convergence map of size N x N into J = jmax + 1 sub-arrays 
of size N x N as follows: 


jmax 


K(x) =cy(xy) + }  wj(x,y), (18) 
j=l 


where jmax represents the number of wavelet bands (or, equivalently, aperture mass maps) 
considered, cz represents a smooth (or continuum) version of the original image x, and wj 


represents the input map filtered at scale 2/ (i.e. the aperture mass map at 0 = 2j pixels). 


Using the wavelet formalism to derive the aperture mass statistic presents a number of 
advantages. Many families of wavelet functions have been studied in the statistical literature, 
and all these wavelet functions could be applied to the aperture mass statistic. This allows us 
to tune our filter function to optimise the signal-to-noise in the resulting maps. In addition, 
for some specific wavelet functions, discrete and very fast algorithms exist, allowing us to 
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compute a set of wavelet scales through the use of a filter bank with a very limited number 
of operations. See [6] for a full review of the different wavelet transform algorithms. 


In the starlet transform algorithm, the wavelet (x,y) is separable and can be defined by: 
* Y) ea pane 
0 (34) =4 [0G TIE (19) 


where (x,y) = ¢(x)g(y) and g(x) is a scaling function from which the wavelet is generated. 
In the case of the starlet wavelet, p(x) is a B3-spline: 


1 
glx) = zy (lx — 27 — 4x — 1° + 6 |x|? — 4]x + 11? + |x + 2°), (20) 


which is a compact function that is identically zero for |x| > 2. 


This wavelet function has a compact support in real space, is well localized in Fourier 
domain, and the wavelet decomposition of an image can be obtained with a very fast 
algorithm (see [24] for a full description). 


Figure 2 shows the starlet wavelet and aperture mass filters defined above in both real and 
Fourier space. Notice that the two filters presented have very similar shapes in real space, 
but different widths. While the starlet filter function goes to zero identically at 0 = 26, 
and remains zero beyond this value, the Map filter function goes to zero as 8 — œ, and 
must therefore be truncated when applied in practice. Clearly at 0 = 50, the Map filters are 
sufficiently close to zero, so that this is an appropriate truncation radius; however this will 
impact the computation time of the Map statistic. 


In Fourier space, we show the response of the Map filter for truncation radii of Oeut = 0 
and 50. Notice that the shape of the response curves in Fourier space for both the Map and 
wavelet filter functions are similar, but that the peak of the response curve for the Map filter 
function shifts when Veut is varied. This is because the effective scale of the filter function 
is being changed. Also notice that in the case of Veut = 0, high-frequency oscillations are 
present in the response curve. This is a direct result of the truncation of the filter function, 
and will occur whenever such truncation is applied. Indeed, imperceptibly small oscillations 
are still present in the response curve for the filter truncated at Veut = 50. 


Such oscillations are not present at all with the starlet filter function, as no truncation is 
applied to this function whatsoever. This gives the starlet transform the distinct advantage 
of being directly applicable, without any consideration needed regarding truncation radii 
and the associated impact on the Fourier-space response of the filter functions. 


3.3. Advantages of the wavelet formalism 


The aperture mass formalism tends to be the preferred method for weak lensing studies for a 
number of reasons. Firstly, the filter functions can be simply expressed analytically, and any 
associated statistics of the aperture mass can therefore be straightforwardly computed from 
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(a) Map and wavelet filter functions in real space (b) Filter response curves in Fourier space 


Figure 2. A comparison of the Map and starlet wavelet filters in real and Fourier space. 


the shear catalogue directly. This avoids the need to generate an aperture mass map, which 


can be computationally intensive, and furthermore does not require the computation of the 
convergence, K. 


Computing the convergence from the shear measurements can be tricky for a number of 
reasons. Equation (4) describes a convolution over all 2D space. Given that we aim to invert 
this equation for an image of a finite size, a direct inversion in Fourier space will give rise 
to significant edge effects, and a leakage of power into so-called B-modes, which effectively 
imply a spurious (non-lensing) cross-component to the shear field y, and usually only arise 
due to systematic effects in the lensing measurements. 


Several methods have been devised to invert equation (4) whilst minimising these 
undesirable effects [26-29]. Most recently, the authors in [30] have presented a method 
based on a wavelet-Helmholtz decomposition, with which they demonstrate a reconstruction 
error at the few percent level, as compared to an error of ~ 30% seen with Fourier-based 
methods. This implies that there is little advantage to the shear catalogue as opposed to 


the convergence map. Indeed, working on the convergence map and under the wavelet 
framework may offer some distinct advantages. 


The first advantage comes directly from the nature of the wavelet transform. The aperture 
mass filter is a bandpass filter. Because the Map reconstruction is only computed for a 
discrete set of aperture scales, information on intermediate scales may be lost, particularly 
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at both the small and large frequency extremes. In contrast, the wavelet transform retains 
information on all scales. The first wavelet scale is effectively a high-pass filter, retaining 
all the high-frequency information in the image, while the remaining wavelet scales are 
bandpass filtered as with the M,p-filtered images. Finally, the wavelet transform retains cy, 
which encodes the large-scale information in the image, and therefore consists of a smoothed 
version of the image (or an image with a low-pass filter applied). Figure 3(a) demonstrates 
this point for a wavelet transform with J = 5, showing the Fourier-space response of the 
wavelet filter as a function of scale j. 


Proces 


Aperture Radius (pixels) 


(a) The filter response of each of 5 scales of the wavelet (b) Comparison of a the processing time of a 
transform brute-force aperture mass algorithm and the starlet 
transform algorithm 


Figure 3. Illustration of the various advantages offered by the wavelet transform in terms of (a) retention of information on all 
scales and (b) algorithm processing time. 


In addition, from a mapping perspective, the starlet wavelet transform algorithm is 
substantially faster than a brute-force computation of the aperture mass statistic in real space. 
Such a naive implementation has complexity « O(N782,,), where N x N is the dimension 
of the image, and Veut is the chosen truncation radius. This scaling means that for large 
apertures or, equivalently, for high-resolution images, the algorithm may prove to be very 
time-consuming. The starlet wavelet transform algorithm is of complexity « O(N?J), where 
J is the number of scales considered, and is limited by N > 2]. This means that the processing 
time for the wavelet transform algorithm is sensitive only to the number of scales considered, 
rather than the size of the filter functions involved, and depends linearly on this number. 


In Figure 3(b), we compare the processing time for the aperture mass algorithm and the 
starlet transform algorithm, both programmed in C++, to analyse an image of 1024 x 1024 
pixels on a 2 x 2.66GHz Intel Xeon Dual-Core processor. We consider aperture scales 0 = 
(2,4, 8,16, 32,64], which correspond to J = jmax +1 = [2,3,4,5,6,7] wavelet scales in the 
wavelet transform. In the aperture mass algorithm, the filters are truncated at a radius 
of cut = 9. For the application of filters which necessitate truncation at a much larger 
radius, we expect the computation time to be roughly an order of magnitude longer. Even 
at the smallest aperture radius, the wavelet transform is ~ 5x faster than the aperture mass 
algorithm. At 0 = 64 pixels, the wavelet transform is ~ 1200x faster than the aperture mass 
algorithm. Note that the wavelet transform for a given value of J simultaneously computes 
the wavelet transform at all scales 2/, 0 < j < J — 1, in addition to the smoothed continuum 
map cy, whilst the aperture mass algorithm computes the transform at a single scale 0. We 
note further that the computational time for the wavelet transform for J = 7 wavelet scales 
is still a factor of ~ 2x less than the computational time for the aperture mass algorithm at 
0 = 2 pixels. 
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As discussed in [31], this time advantage further extends to higher-order statistics of the 
Map, which are typically related directly to n-point correlation functions of the shear. In 
recent years, tree codes have been employed to speed up computation of n-point correlation 
functions. Typical tree codes to compute n-point correlation functions are O(Ngqi log(Neal)) 
[32] on a shear catalogue. For a single Euclid exposure of 0.5 deg’, we can expect Neal ~ 
54,000 (30 galaxies/ arcmin?) - 180,000 (100 galaxies / arcmin2). Tree codes exist that act on 
pixelated data [33] which run at O (Noix!bin) where Nein is the total number of pixels and 
Npin is the number of bins in the correlation function. For a Euclid exposure, assuming pixels 
of 1 arcminute, we have N2.. = 1800, and npin will be dependent on the required resolution 


pix 
of the correlation function. 


The wavelet method acts on pixellated data, and is O(N? J) in computation time, so our 
algorithm will be comparable for computation of 2-point statistics, if the 2-point correlation 
function is computed on pixellated data, but if a shear catalogue is used, we will have a 
faster algorithm by at least an order of magnitude. For higher-order statistics, this advantage 
is even more pronounced. Furthermore, while optimised software is freely and publicly 
available to compute the wavelet transform, such optimised software is not available for 
n-point correlation functions. 


Another advantage of the wavelet formalism is that it is possible to carry out an explicit 
denoising of the convergence field using thresholding based on a False Discovery Rate 
method. For details on this method, see [34], where the MRLens software package encoding 
this method is presented. This allows one to derive robust detection levels in wavelet space, 
and to produce high-fidelity denoised mass maps. 


In addition, wavelet-based methods offer more flexibility than aperture mass filters. Whilst 
we have thus far discussed only the starlet wavelet function, many other wavelet dictionaries 
may be used. The starlet filter seems ideal for lensing studies, due to its similarity to the 
aperture mass filter presented here, which was deemed to be optimal in [22]. However, 
different dictionaries may be optimal in different applications; for example, if one were 
attempting to study filamentary structure, ridgelets or curvelets might be a more appropriate 
basis. The vastness of the wavelet libraries and the public availability of fast algorithms to 
compute these transforms are major strengths of wavelet-based approaches. 


4. Weak lensing beyond two-point shear statistics 


As noted previously, the optimality of second-order statistics to constrain cosmological 
parameters depends heavily on the assumption of Gaussianity of the field. However the 
weak lensing field is composed, at small scales, of non-Gaussian features such as clusters 
of galaxies. These non-Gaussian signatures carry additional information that cannot be 
extracted with second-order statistics. Many studies [35-39] have shown that combining 
second-order statistics with higher-order statistics tightens the constraints on cosmological 
parameters. We now consider the higher-order statistics most commonly used in weak 
lensing studies aimed at detecting and constraining non-Gaussianity in the lensing field. 


3 The MRLens software, along with many other software packages for astronomical applications, is freely available 
here: http:/ /cosmostat.org/software.html 
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4.1. Higher-order lensing statistics 


As we have already noted, it is more convenient to consider higher-order statistics of the 
Map or wavelet transform of the convergence field, as opposed to statistics of the shear field 
directly. The obvious first extension to two-point statistics is to consider the three-point 
correlation function or, equivalently, the bispectrum, usually considered as a function of 
aperture or wavelet scale. 


Interesting analytical results relative to the shear three-point correlation function or the 
convergence bispectrum have been reported (for example, see [40-43]). However, it has 
been shown that in using only the equilateral configuration of the bispectrum, the ability to 
discriminate between cosmological models is relatively poor [38]. An analytical comparison 
has been performed in [39] between the full bispectrum and an optimal match-filter peak 
count for a Euclid-like survey, and both approaches were found to provide similar results. 
However, as the full bispectrum calculation has a much higher complexity than other 
statistics, and no public software exists to compute it, we will not consider the full bispectrum 
here. Rather, we will restrict ourselves to statistics that are more straightforward to compute 
from the Map or wavelet transform maps. 


4.1.1. The Skewness 


The skewness of the aperture mass map, (Min), is the third-order moment of the aperture 
mass Map (0) and can be computed directly from shear maps filtered with different aperture 
mass. The skewness is a measure of the asymmetry of the probability distribution function. 
The probability distribution function will be more or less skewed positively depending on the 
abundance of dark matter haloes at the @ scale. The formalism exists to predict the skewness 
of the aperture mass map for a given cosmological model, which is related to the three-point 
correlation function or the bispectrum of density fluctuations 6 [20, 37]. 


In [37], it is argued that the skewness as a function of scale is a preferable statistic to 
three-point correlation functions of the shear, as the integral relations between (Mi,) and the 
bispectrum are much easier and faster to compute than the three-point correlation function. 
This is because the skewness is a local measure of the bispectrum, whereas the integral kernel 
for the three-point correlation function is a highly oscillating function with infinite support. 


This statistic can equivalently be computed from the wavelet transform of the convergence 
map as the skewness (w? ) of the wavelet band j corresponding to the aperture scale 0. This 
can be computed either on the noisy convergence map estimated, for example, by the Fourier 
space relationship between the shear and convergence — in which case the results should be 
comparable to the Map results — or on a denoised convergence map using the method of [34], 
which should provide improved results. 


4,1.2. The Kurtosis 


The kurtosis of the aperture mass map, (Map), is the fourth-order moment of the aperture 
mass Map (8) and can be computed directly from the different aperture mass maps. The 
kurtosis is a measure of the peakedness of the probability distribution function. The presence 
of dark matter haloes at a given 0 scale will flatten the probability distribution function and 
widen its shoulders leading to a larger kurtosis. The formalism exists to predict the kurtosis 
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of the aperture mass map for a given cosmological model, which is related to the four-point 
correlation function or trispectrum of the 3D density contrast [20]. 


Again, this statistic can be computed from the wavelet transform of the convergence map 
as the kurtosis (w4) of the wavelet band j corresponding to the aperture scale 0. This can 
be computed either on a noisy or a denoised convergence map, the latter expected to yield 
improved constraints. 


4,1.3. Peak Counts 


We define a peak as set of connected pixels above a detection threshold 7, and denote the 
peak count in the Map and wavelet maps by Ph, and P respectively. If peak counting 
is carried out on a denoised map, the detection threshold 7 is automatically set by the 
denoising algorithm, and a small threshold € is used when identifying peaks in the denoised 
maps, in order to reject spurious detections in these denoised maps. 


We consider all pixels that are connected via the sides or the corners of a pixel as one 
structure. In a two-dimensional projected map of the convergence, we are therefore unable to 
discriminate between peaks due to massive halos and peaks due to projections of large-scale 
structures such as filaments. 


While theory exists to predict cluster counts from the halo model and a cosmological 
model encoding the growth and evolution of structures in the Universe, there is no analytic 
formalism to predict the fraction of spurious detections in lensing maps arising from 
projections of large-scale structure. In [44], an attempt is made to derive an analytic 
formalism for predicting peak counts in projected maps. However, this method is based 
on an assumption of Gaussianity in the lensing field and, predictably, underestimates counts 
at the high end of the mass function, where clusters of galaxies are dominant. This means 
that predicting peak counts for a given cosmological model relies on considering N-body 
simulations generated under a range of cosmological parameters (for example, see [45]). 


4.2. Optimal capture of non-Gaussianity 


The question now arises: which of these statistics provides the most information about 
the underlying cosmology and - specifically — non-Gaussianity in the density field? In 
order to asses the performance of these statistics, we consider the ability of each statistic 
to discriminate between cosmological models using N-body simulations under a range of 
different cosmologies [38, 46]. 


To this end, we consider N-body simulations carried out for 5 different cosmological models 
along the Om — og degeneracy, the parameters of which are summarised in Table 1 below 
[46]. These simulations were carried out using the RAMSES N-body code [47], and full details 
of the models considered are given in [38]. For each model, 100 realizations of weak lensing 
maps were generated for a field of 3.95° x 3.95°, downsampled to 1024 x 1024 pixels (0.23/ 
per pixel). Noise was added to the simulations at a level consistent with predictions for deep 
space-based observations, with a number density of galaxies of n g = 100 gal/ arcmin?. This 
is somewhat optimistic; however, as seen in [46], increasing the noise level in the simulations 
does not change the conclusions of the comparison test between the higher-order statistics. 
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Model Box Length Om On h Tg 
(h- Mpc) 
mı 165.8 0.23 0.77 0.594 1 
mz 159.5 0.3 0.7 0.7 0.9 
m3 152.8 0.36 0.64 0.798 0.8 
m4 145.7 0.47 0.53 0.894 0.7 
ms 137.5 0.64 0.36 0.982 0.6 


Table 1. Parameters of the five cosmological models that have been chosen along the Om — og degeneracy. The simulations 
have 256° particles, h is equal to Hy /100kms~! Mpc~?, where Hp is the Hubble constant. 


To find the best statistic, we need to characterise quantitatively for each statistic the 
discrimination between two different models mı and mọ. To do this, we consider the 
distribution functions for the different statistics estimated on the 100 realisations of each 
model. These distributions are expected to overlap, so False Discovery Rate (FDR) method 
is used to determine the threshold 7, in the distribution function of model m4, such that the 
distribution function of model m) accounts for fewer than a fraction æ = 0.05 of the total 
counts. A similar threshold t is defined for the distribution function of model m2. This 
is illustrated in Figure 4. The discrimination efficiency is then defined as the percentage of 
counts remaining under the hatched area for each model. 


Figure 4. Illustration of the discrimination efficiency criterion. t} and T2 are defined such that the yellow area under the curve 
represents 0.05 of the total area under the curves for models mz and mı, respectively. The discrimination efficiency is defined 
as the area under the m; curve delimited by the threshold 7;, as a percentage of the total area under the model m; curve. 


In Table 2 below, we show the mean discrimination efficiency obtained considering all pairs 
of models for the higher-order statistics of interest. These are presented as a function of 
Map or wavelet scale, and shown for statistics computed on the aperture mass map and the 
wavelet transform of the noisy convergence map. The results overall are comparable, with the 
wavelet transform statistics appearing to offer a slight improvement over the Map statistics 
in all cases. This is perhaps related to the fact that, when computing the Map statistics, we 
truncate the filter at cut = 8, which may give rise to some small systematics in the resulting 
maps. Furthermore, it is clear that, in all cases, the peak statistics are much more efficient at 
discriminating between cosmological models than either the skewness or the kurtosis. 
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“pa. (nfo\ (Mt) pe pao. 
0; (Map) (Map) PM, Př, 


0.46 04.60% 0230% 39.70% 54.95 % 
0.92 33.40% 03.70% 79.30% 76.40% 
1.85’ 03.45% 01.45% |91.25%| 89.20 % 
3.70 15.15% 23.00% 69.40% 86.70% 
740 26.95% 24.30% 490% 60.50% 


(a) Mean discrimination efficiency of Map statistics 


Scale (ws ) (w ) P, 0 L w; 


J J 
046 0200% 01.15% 12.05% 00.70% 
0.92 37.95% 0475% 8630% 73.05% 
1.85’ 03.55% 02.10% [94.40%] 93.85 % 
3.70 18.25% 25.65% 8405% 87.05% 
740 36.40% 30.90% 24.60% 66.35% 


(b) Mean discrimination efficiency of wavelet statistics 


Table 2. Mean discrimination efficiency (in percent) from noisy Map reconstructions of the shear field, compared with that 
obtained on the wavelet transform of the noisy convergence maps. 


In Table 3 below, we present the mean discrimination efficiency of statistics computed 
on the denoised convergence maps, as well as the peak discrimination efficiency between 
each pair of models in our sample. There is a clear improvement on the discrimination 
efficiency of all the higher-order statistics when denoising is applied to the convergence 
maps. This is unsurprising, as the presence of Gaussian (or near-Gaussian) noise within the 
shear and convergence maps will make the whole field appear more Gaussian, masking the 
non-Gaussian features in the map and pushing the skewness and kurtosis values closer to 
zero. The more noise that is present in the maps, the stronger this effect will be. It is for 
this reason that such a dramatic improvement is seen in these statistics when denoising is 
applied to the convergence maps before the statistics are computed. 


The improvement seen in the discrimination efficiency of the peak statistics is also significant, 
and this arises due to the fact that while the peak counts in the noisy maps involve application 
of a simple no detection threshold to the maps, the denoising algorithm involves a more 
sophisticated discriminant between signal and noise, again making use of the FDR method. 
This technique is more effective and distinguishing between signal and noise, and therefore 
more information about the true peaks in the map is retained using this denoising method 
as compared to the no thresholding applied previously. 


It is clear that peak statistics are much more efficient at capturing non-Gaussianity in the 
shear field than other higher-order statistics, and that it is advantageous to measure this 
statistic on denoised convergence maps. However, in order to combine constraints from 
peak counting with other probes, such as two-point statistics, in order to better constrain 
the cosmological model, it is important to be able to predict peak statistics as a function 
of cosmology in order to be able to compute the likelihood function for all models within 
the parameter space. As we have discussed, this is not possible to do analytically, due to 
contamination by projections of large scale structures, and it is too computationally-intensive 
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Scale (@) (a) Po, 


0.46 53.40% 43.20% 68.35% 
0.92 47.90% 4115% 9245 % 
1.85 58.80% 44.70% |96.75 % 
3.70 63.30% 48.05% 90.40% 
740 5490% 4045% 6345% 


(a) Mean discrimination efficiency 


mı m2 m3 Ma ms5 
mı x 85% 100% 100% 100% 
mı 89% x 92% 100% 100% 
m3 100 % 92 % X 89 % 100 % 
mą 100 % 100 % 92 % x 98 % 
ms 100% 100% 100% 98% x 


(b) Discrimination efficiency for all models obtained using peak statistics 


Table 3. Discrimination efficiency for statistics computed on denoised convergence maps. 


to consider obtaining the predictions, and associated covariances, from N-body simulations 
if we wish to sample the parameter space completely and at high resolution. 


As the quality of data available to astronomers continues to improve, it has recently become 
possible to consider the shear field in three dimensions, using the colours of galaxies to 
determine their redshifts and, therefore, distances from us. If this information may be used 
to deproject the lensing signal, and thus recover information about the full three-dimensional 
density field, it should then be possible offer predictions for peak statistics as a function of 
cosmology, uncontaminated by projection effects which will, in turn, allow us to improve our 
constraints on our cosmological model. 


5. Reconstructing the density contrast in 3D 


The measured shear can be related to the 2D convergence via a simple linear relationship 
(equation (4)), and the convergence is related to the density contrast by another linear 
mapping (equation (6)). We can express these relationships conveniently in matrix notation 
as 


a 
— 
© 
N 
ys 
II 
ie) 
S 
= 
Ss 
N 


(21) 
y(0,z) = PyxQ6(0,z) . (22) 


The 3D lensing problem is therefore one of finding an estimator to invert equation (21) 
or (22) in the presence of noise. For simplicity, in what follows we assume the shear 
measurement noise is Gaussian and uncorrelated between redshift bins. In practice, errors 
in the photometric estimates of the redshifts of galaxies will introduce correlations between 
redshift bins. Methods have been developed to account for such errors, however, and the 
problem is therefore readily tractable [48]. 
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Therefore the 3D lensing problem is effectively one of observing the density contrast 
convolved with the linear operator R, and contaminated by noise, which is assumed to be 
Gaussian. Formally, we can write 


d=Rs+e, e~N(0,07), (23) 
where d is the observation, s the real density and € the Gaussian noise. 


5.1. Linear approaches to 3D map-making 


The general idea behind linear inversion methods is to find a linear operator H that acts on 
the data vector to yield a solution that minimises some functional, such as the variance of 
the residual between the estimated signal and the true signal, subject to some regularisation 
or prior-based constraints. Two different linear approaches have been described in recent 
literature [4, 49]. 


The most competitive method is that proposed in [4], in which the authors propose a 
Saskatoon filter [50, 51], which combines a Wiener filter and an inverse variance filter, with 
a tuning parameter « introduced that allows switching between the two. This gives rise to a 
minimum variance filter, expressed as 


éuv = [41+ SRE İR] SR'S d, (24) 


where S = (sst) encodes prior information about the signal covariance, & = (nn*) gives 
the covariance matrix of the noise, and 1 is the identity matrix. 


This switching is designed to allow a balance between the increased constraining power 
offered by the Wiener filter over the inverse variance filter - which yields an improved 
signal-to-noise in the reconstruction — and the biasing that the Wiener filter imposes on the 
solution. 


As discussed extensively in [4] and [49], linear methods give rise to a significant bias in the 
location of detected peaks, damping of the peak signal, and a substantial smearing of the 
density along the line of sight. Furthermore, the resolution attainable in the reconstructions 
obtained using linear methods is limited to the resolution of the data. In other words, we 
cannot reconstruct the density contrast field at higher resolution than the resolution of our 
input data which, in turn, is limited by the noise properties of the data. 


Figure 5 shows the 3D reconstruction obtained in this way for a simulated cluster of 
galaxies at a redshift of zą = 0.25. The tuning parameter used for this reconstruction was 
« = 0.05. The cluster was simulated according to an NFW halo with M = 10471 Mo 
and c = 3, and the shear data were assumed to come from a galaxy distribution given by 


p(z) « ze (14z)!° [52, 53], with a maximum redshift of Zmay = 2.0 and a galaxy density 
of ng = 100 galaxies/ arcmin?. The simulation covers a 1° x 1° field binned into 60 x 60 
angular pixels, and 20 redshift bins. Shown in the figure are both the 3D rendering of the 
reconstruction, and a 1D plot showing the four central lines of sight through the cluster as a 
function of redshift. The smearing, damping and redshift bias effects are all clearly visible in 


the 1D plot, where the amplitude of the cluster density contrast should be ô ~ 36. 
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Radial Wiener Filter 
T T 
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(a) 3D plot of the cluster reconstruction. (b) The four central lines of sight 


through the cluster. 


Figure 5. The reconstruction of a simulated cluster of galaxies using the method of [4]. Coloured broken lines in the 1D plot 
show the reconstruction of the four central lines of sight, while the solid line shows the true cluster density contrast. 


5.2. Compressive sensing 


We consider some data Y; (i € [1,..,m]) acquired through the linear system 
Y=0X, (25) 


where © is an m x n matrix. Compressed sensing [54, 55] is a sampling/compression theory 
based on the sparsity of the observed signal, which shows that, under certain conditions, one 
can exactly recover a k-sparse signal (a signal for which only k pixels have values different 
from zero, out of n total pixels, where k < n) from m < n measurements. 


This recovery is possible from undersampled data only if the sensing matrix © verifies 
the restricted isometry property (RIP) [see 54, for more details]. This property has the effect 
that each measurement Y; contains some information about the entire pixels of X; in other 
words, the sensing operator © acts to spread the information contained in X across many 
measurements Yj. 


Under these two constraints — sparsity and a transformation meeting the RIP criterion — a 
signal can be recovered exactly even if the number of measurements m is much smaller than 
the number of unknown n. This means that, using CS methods, we will be able to outperform 
the well-known Shannon sampling criterion by far. 


The solution X of (25) is obtained by minimizing 


min |X|], sf Y=0X, (26) 


where the £; norm is defined by ||X||; = ); | X; |. The 4 norm is well-known to be a 
sparsity-promoting function; i.e. minimisation of the /4 norm yields the most sparse solution 
to the inverse problem. Many optimisation methods have been proposed in recent years to 
minimise this equation. More details about CS and ¢; minimisation algorithms can be found 
in [6]. 
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In real life, signals are generally not “strictly” sparse, but are compressible; i.e. we can 
represent the signal in a basis or frame (Fourier, Wavelets, Curvelets, etc.) in which the curve 
obtained by plotting the obtained coefficients, sorted by their decreasing absolute values, 
exhibits a polynomial decay. Note that most natural signals and images are compressible in 
an appropriate basis. 


We can therefore reformulate the CS equation above (Equation (26)) to include the data 
transformation matrix ®: 


min lali st. Y=O@a, (27) 


where X = ®*q, and a are the coefficients of the transformed solution X in ®, which 
is generally referred to as the dictionary. Each column represents a vector (also called an 
atom), which ideally should be chosen to match the features contained in X. If ® admits a 
fast implicit transform (e.g. Fourier transform, Wavelet transform), fast algorithms exist to 
minimise Equation (27). 

One problem we face when considering CS in a given application is that very few matrices 
meet the RIP criterion. However, it has been shown that accurate recovery can be obtained 
, is low [56]. 
The mutual coherence measures the degree of similarity between the sparsifying basis and 


the sensing operator. Hence, in its relaxed definition, we consider a linear inverse problem 
Y = O®X as being an instance of CS when 


as long the mutual coherence between © and ®, pọ, = max; (0; ®,, ) 


1. the problem is underdetermined, 
2. the signal is compressible in a given dictionary ®, 


3. the mutual coherence pig@,o is low. This will happen every time the matrix A = @® has 
the effect of spreading out the coefficients «j of the sparse signal on all measurements Yj. 


Most CS applications described in the literature are based on such a soft CS definition. 
Compressed sensing was introduced for the first time in astronomy for data compression 
[57, 58], and a direct link between CS and radio-interferometric image reconstruction was 
recently established in [59], leading to dramatic improvement thanks to the sparse 44 recovery 
[60]. 


The 3D weak lensing reconstruction problem can be seen to completely meet the soft-CS 
criteria above. The problem is underdetermined, as we seek a higher resolution than can 
be attained in the noise-limited observations, the matter density in the Universe is sparsely 
distributed, and the lensing operator spreads out the underlying density in a compressed 
sensing way. 


In particular, for the reconstruction of clusters of galaxies, we are in a perfect situation for 
sparse recovery because clusters are localised in the angular domain, and are not resolved 
along the line of sight owing to the bin size. They can therefore be modelled as Dirac 
6—functions along the line of sight, while an isotropic wavelet basis can be used in the 
angular domain. 
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5.3. Results and future prospects 


In [61], we present an algorithm to solve the 3D lensing problem. In this method, the 3D 
lensing problem is reduced to a one-dimensional problem, by taking as the data vector the 
(noisy) lensing convergence along each line of sight, which is related to the density contrast 
through Equation (21). Each line of sight can therefore be considered independently. 


Clearly, a one-dimensional implementation throws away information, because we do not 
account at all for the correlation between neighbouring lines of sight that will arise in 
the presence of a large structure in the image; however, reducing the problem to a single 
dimension is fast and easy to implement, and allows us to test the efficacy of the algorithm 
using a particularly simple basis function through which we impose sparsity. 


In Figure 6, we show the reconstruction obtained in this way for the simulated cluster 
described in section 5.1. The line of sight plot shows a clear improvement in all the 
target areas, reducing the bias, smearing and damping effects seen using linear methods. 
Small-scale noise is present, particularly at high redshifts, but this is likely due to overfitting 
in the algorithm used, and may be reduced by considering the whole 3D field, rather than 
each line of sight independently (for a full discussion of this issue, see [61]). 


(a) 3D plot of the cluster reconstruction. (b) The four central lines of sight 
through the cluster. 


Figure 6. The reconstruction of a simulated cluster of galaxies using the method of [61]. Coloured broken lines in the 1D plot 
show the reconstruction of the four central lines of sight, while the solid line shows the true cluster density contrast. 


This reconstruction was undertaken using the same resolution on the reconstruction as on the 
input data. However, we know that the CS approach is particularly well-adapted to dealing 
with ill-posed inversion problems. In order to test this, we consider a cluster at a redshifts 
of Zą = 0.2, simulated as before. We use Nsp = 20 redshift bins in our data, but now aim to 
reconstruct our density contrast with a redshift resolution of Nip = 25. The results are shown 
in Figure 7. Here, again, we see small scale noise at high-redshift, but the overall smearing 
and redshift bias issues seen in the linear reconstructions is absent. 


As noted previously, the one-dimensional solver employed here throws away a wealth 
of information about the angular correlation of the lensing signal, and is thus not 
optimal. Indeed, a simple algorithm based on this CS approach, but implemented as 
a full three-dimensional treatment, offers marked improvements in the quality of the 
reconstructions [62]. 


Exploring the Components of the Universe Through Higher-Order Weak Lensing Statistics 
http://dx.doi.org/10.5772/51871 


Redshift 


(a) 3D plot of the cluster reconstruction. (b) The four central lines of sight 
through the cluster. 


Figure 7. The reconstruction of a simulated cluster of galaxies at a redshift of za = 0.2, and with improved resolution on the 
reconstruction as compared to the input data. 


This is demonstrated in Figure 8, which shows reconstructions of the cluster at redshift 
Zq = 0.2, with the same improved resolution in the reconstruction as before, but this time 
using the three-dimensional CS approach. Dramatic improvement is seen on all fronts, with 
the reconstructions showing no bias or redshift smearing, and very little amplitude damping, 
and with none of the small-scale false detections seen in the 1D CS approach. 


r 
7 aaa 
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(a) 3D plot of the cluster reconstruction. (b) The four central lines of sight 
through the cluster. 


Figure 8. The reconstruction of a simulated cluster of galaxies at a redshift of za = 0.2 using a 3D CS approach, and with 
improved resolution on the reconstruction as compared to the input data. 


This marked improvement in the cluster reconstruction seen in the 3D CS approach 
represents a definite step in the right direction for weak lensing studies. There is much 
work to be done, however. The application of this CS approach to the 3D lensing problem is 
a very recent development, and many questions remain: the choice of algorithm, for example; 
how best to control the noise in the reconstruction; how to deal with photometric redshift 
errors; by what factor the reconstruction resolution might be improved as compared to the 
data. 
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Yet it is clear that this approach opens up the possibility of being able to generate accurate 
reconstructions of the density contrast using weak lensing measurements, and perhaps using 
information such as the 3D peak count — in combination with constraints from other probes — 
to place ever-tighter constraints on our cosmological model. This, in turn, will offer a unique 
insight into the nature and properties of the dark components of the Universe. 
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